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ABSTRACT 

Cosmological inference with large galaxy surveys requires theoretical models that combine precise predictions for large-scale 
structure with robust and flexible galaxy formation modelling throughout a sufficiently large cosmic volume. Here, we introduce 
the MILLENNIUMTNG (MTNG) project which combines the hydrodynamical galaxy formation model of ILLUstrisTNG with 
the large volume of the MiLLENNIUM simulation. Our largest hydrodynamic simulation, covering (500 h~'Mpc)? ~ (740 Mpc)’, 
is complemented by a suite of dark-matter-only simulations with up to 4320° dark matter particles (a mass resolution of 
1.32 x 108 h-'Mo) using the fixed-and-paired technique to reduce large-scale cosmic variance. The hydro simulation adds 43207 
gas cells, achieving a baryonic mass resolution of 2 x 10’ h~'Mo. High time-resolution merger trees and direct lightcone outputs 
facilitate the construction of a new generation of semi-analytic galaxy formation models that can be calibrated against both 
the hydro simulation and observation, and then applied to even larger volumes — MTNG includes a flagship simulation with 
1.1 trillion dark matter particles and massive neutrinos in a volume of (3000 Mpc)?. In this introductory analysis we carry out 
convergence tests on basic measures of non-linear clustering such as the matter power spectrum, the halo mass function and halo 
clustering, and we compare simulation predictions to those from current cosmological emulators. We also use our simulations 
to study matter and halo statistics, such as halo bias and clustering at the baryonic acoustic oscillation scale. Finally we measure 
the impact of baryonic physics on the matter and halo distributions. 
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1 INTRODUCTION sions such as Euclid (Laureijs et al. 2011) and Roman (Spergel et al. 
2015), as well as new powerful telescopes on Earth such as DESI 
(DESI Collaboration 2016), PFS (Takada et al. 2014) and Rubin 
(LSST Science Collaboration 2009), will map out billions of galax- 
ies through extremely large regions of space. They will primarily use 
various measures of galaxy clustering and weak gravitational lensing 
to carry out meticulous tests of the cosmological model. The primary 
goals are to search for deviations of dark energy from a cosmological 
constant, for non-gaussianities in the primordial fluctuation field, for 
signatures of a law of gravity different from general relativity, and to 
measure the mass of the light neutrino flavors. 


The amazing progress in observational cosmology over the last 
decades has brought many surprises. Perhaps the most stunning is 
that we live in a Universe where most of the matter is comprised of 
yet unidentified collisionless dark matter particles, while ordinary 
baryons produced in the Big Bang make up only a subdominant 
part (Aghanim et al. 2020). Furthermore, in the last 5 billion years 
or so, a dark energy component has progressively become stronger 
and begun to overwhelm the matter density, driving an accelerated 
expansion of the Universe. The real physical nature of dark energy, 
the identity of the dark matter particles, as well as the mass of the 


neutrinos which contribute a tiny admixture of hot dark matter, are 
profound and fundamental open questions in physics. 

To make further progress, the firmly established ACDM standard 
cosmological model will be subjected to precision tests in the coming 
years that are far more sensitive than anything done thus far. Forth- 
coming cosmological mega galaxy surveys carried out by space mis- 
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It has widely been recognised that systematic uncertainties in our 
ability to compute very accurate theoretical predictions for ACDM, 
as well as neighbouring cosmological models, could become a limit- 
ing factor in making full use of the statistical power of the upcoming 
data. Simulation predictions need to cover much larger cosmologi- 
cal volumes than customary so far to match the statistical power of 
the new surveys. They also have to be able to predict the non-linear 
matter clustering highly accurately, ideally account for the impact 
of baryonic physics on matter clustering in a reliable fashion, and 
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produce realistic galaxy properties and galaxy clustering signals. 
In addition, one would like to be able to create such predictions 
not only for one set of cosmological parameters, but for many model 
variants in a computationally inexpensive fashion. Reaching this goal 
is very demanding, and will likely require an innovative combina- 
tion of “ground-truth” simulations based on full-hydrodynamical and 
N-body simulations (see Vogelsberger et al. 2020, for a review), ap- 
proximate but fast simulation techniques (e.g. Feng et al. 2016), 
rescaling techniques (Angulo & White 2010), semi-analytic galaxy 
formation methods (e.g. White & Frenk 1991, SAM), more empir- 
ical methods such as halo occupation distribution modelling (e.g. 
Berlind & Weinberg 2002, HOD) and subhalo abundance match- 
ing (e.g. Conroy et al. 2006, SHAM), plus data-driven approaches 
such as machine-learning techniques (e.g. Villaescusa-Navarro et al. 
2021). 

A number of groups have in recent years produced very large sim- 
ulation models as initial steps to tackle this problem. Among them 
are the MILLENNIUM-XXL (Angulo et al. 2012), the Euclid FLAGsHip 
(Potter et al. 2017), the OurERRm™ (Heitmann et al. 2019), the ABa- 
cusSumnmit (Maksimova et al. 2021), or the Ucnuu (Ishiyama et al. 
2021) simulations. These simulations are characterised by their large 
cosmological volume (= 23 h~3Gpc3) and by their large number of 
resolution elements (Vp 2 6000), making them ideally suited to 
match the statistical power of the large-scale galaxy surveys. How- 
ever, most of these simulations do not employ a physically motivated 
galaxy formation model to produce their mock catalogues. Instead, 
those mocks are often created by empirical models, such as HOD, 
which do not take into account the impact of baryonic physics on the 
internal structure of haloes and the clustering of matter. Moreover, 
most of these large simulation projects do not provide detailed dark 
matter (sub)halo merger trees which are needed to create realistic 
mocks through semi-analytical models of galaxy formation. 

In this paper, we introduce a new project along this line of research, 
which we have named “MiLLENNIUMTNG” based on its close connec- 
tions to two older simulation projects, the MILLENNIUM simulation 
(Springel et al. 2005), and The Next Generation Ilustris Simulations 
(ILLustrisTNG; Pillepich et al. 2018b; Nelson et al. 2018; Springel 
et al. 2018; Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 
2019b; Pillepich et al. 2019). Both have managed to substantially 
advance galaxy formation modelling, Millennium by introducing the 
first 10 billion particle simulation in a 500 h-'Mpc box that was aug- 
mented with subhalo merger trees that allowed the construction of 
detailed semi-analytic galaxy formation models, while IlustrisTNG 
excelled by combining an accurate moving-mesh hydrodynamical 
technique with a sophisticated model for galaxy formation physics, 
as well as the use of a set of different box size, TNG50, TNG100, and 
TNG300. Still, even TNG300 has a boxsize! of only 205 h~!Mpc 
— clearly too small for the required precision on large cosmological 
scales. 

In the MILLENNIUMTNG project we push the hydrodynamical 
modelling of TNG to a volume nearly 15 times larger, reaching the 
500 h-'Mpc (= 740 Mpc) ona side that could be done with N-body 
techniques in the Millennium, more than a decade ago. To achieve 
this we use an unprecedentedly large number of resolution elements 
for a high-resolution hydrodynamical simulation of galaxy forma- 
tion, with a mass resolution that is still very close to that of TNG300 
and thus sufficient to model the formation of large galaxies with rea- 
sonable accuracy. We accompany this simulation with a set of dark 


! Which is very close to 300 Mpc without the conventional h~!, hence the 
name TNG300. 
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matter only simulations, in the same volume and with the same initial 
phases. While we use the same volume, our best mass resolution is 
nearly an order of magnitude better than that of the original Millen- 
nium simulation. We have also considerably refined subhalo finding 
and tracking, thereby supporting more accurate semi-analytic mod- 
elling. In addition, we carry out these simulations in pairs using a 
variance suppression technique in order to boost the effective volume 
even further. Finally, we augment our simulation set with additional 
N-body runs that take massive neutrinos explicitly into account. Here 
we push also the volume and particle number further, reaching more 
than a trillion particles, and a volume of (3 Gpc)>. 

Our overarching goal with this simulation set is to better link large- 
scale structure studies with non-linear galaxy formation simulations. 
We can achieve this by comparing the full-hydrodynamical simula- 
tions with semi-analytic models applied to the dark matter-only sim- 
ulations, thereby assessing and improving the modelling uncertainty, 
and then by rolling out the semi-analytic model to larger volumes, 
as realised, for example, in our neutrino simulations. Furthermore, 
the simulation set of MTNG provides many opportunities to test 
the internal consistency of the simulation predictions, in particular 
through detailed convergence tests, tests of box size dependencies, as 
well as the ability to assess the impact of baryonic physics and finite 
neutrino masses on matter and galaxy clustering. 

This paper is one of a set of introductory papers we have prepared 
for the MILLENNIUMTNG project. In the present work, we introduce 
the technical aspects of the simulations and present a high-level anal- 
ysis of the matter and halo statistics. Pakmor et al. (2022) describe in 
detail the full-physics MTNG simulation and give a first impression 
of cluster cosmology with MTNG. Barrera et al. (2022) introduce an 
updated version of the L-GALAxies semi-analytic model of galaxy 
formation (Henriques et al. 2015) and its application to the MTNG 
lightcones. Ferlito et al. (2022) analyse weak-lensing convergence 
maps from both DM-only and full-physics runs, while Bose et al. 
(2022) present a galaxy clustering study based on colour-selected 
(blue and red) galaxy samples. Hadzhiyska et al. (2022a,b) present 
an improved halo occupation model (refining the one-halo and two- 
halo terms) of luminous-red and emission-line galaxies using the 
MTNG simulations, and Delgado et al. (2022) study intrinsic align- 
ments of galaxy shapes and compare predictions between the dark 
matter only and full-physics simulations. Kannan et al. (2022) in- 
vestigate properties of very high redshift galaxies (z > 8) in the 
MTNG full-hydrodynamical run. Finally, Contreras et al. (2022) use 
the MTNG simulations and SAM catalogues to infer cosmological 
parameters of SDSS-like samples. These introductory papers cover 
a range of interesting astrophysical and cosmological topics that can 
be addressed with the simulations, however they are far from exhaus- 
tive. We therefore plan to make the data fully publicly available for 
general use by the scientific community in due time, following the 
successful and productive examples of the MILLENNIUM (Lemson & 
Virgo Consortium 2006) and ILLustrisTNG (Nelson et al. 2019a) 
projects. 

The paper is organised as follows. In Section 2, we describe the 
technical aspects and specifications of the MTNG simulations, while 
in Section 3 we introduce the various data sets produced by the 
different simulations. In Section 4, we present results for the foun- 
dational matter and halo statistics, while in Section 5 we extend this 
to an analysis of large-scale (baryonic acoustic oscillations scale) 
matter and halo clustering, as well as scale-dependent halo bias. In 
Section 6, we briefly discuss the baryonic impact on basic matter 
and halo statistics and make a comparison with the results from the 
ILLustRisTNG simulations. Finally, in Section 7 we summarise the 
results and present our conclusions. 
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2 THE MILLENNIUM-TNG SIMULATIONS 


The MILLENNIUMTNG (MTNG) simulations combine large dark 
matter-only and full-physics hydrodynamical computations with the 
goal to link predictions for the evolution of large-scale structure to 
non-linear galaxy formation, while at the same time offering suf- 
ficient volume to allow accurate cosmological inferences. In this 
section, we give an overview of technical specifications of the simu- 
lation set. 

The dark matter-only simulations of the project were run with a 
slightly customised version of the modern GapcET-4 code (Springel 
et al. 2021). Our code extensions relative to the public version are 
primarily concerned with the inclusion of relativistic matter-energy 
components, such as massive neutrinos and a photon background. 
The hydrodynamical simulations have been carried out with the 
Arepo code (Springel 2010; Pakmor et al. 2016; Weinberger et al. 
2020) instead, which has been augmented by us with the group find- 
ing and lightcone outputting routines of the GapGEt-4 code, and has 
furthermore been substantially modified to yield better memory effi- 
ciency and improved scalability when using a very large number of 
processor cores. 

The dark matter simulations consist of one series of runs which all 
use the same volume, i.e. (500 h~!Mpc)? ~ (740 Mpc)3, but which 
vary the number of DM particles systematically from 270? to 43203, 
spaced by a factor of 8 in mass resolution. The highest-resolution run 
(abbreviated MTNG740-DM in the following) improves the original 
MILLENNIUM mass resolution thus by about a factor of 8, which 
translates to a DM particle mass of 1.32 x 108 A“1Mo. Actually, 
the exact value of the ratio of the particle masses is not precisely 8 
because the cosmological parameters we use have changed compared 
to constraints at the time of the Millennium simulation. 

We employ the variance suppression technique of Angulo & 
Pontzen (2016) and run two realisations of each of the DM setups, 
with mode amplitudes set in each case to the square root of the power 
spectrum, and with the second realisation having phases that are mir- 
rored relative to those of the first one. We refer to these runs as the A- 
and B-series of otherwise identical simulations. With such a pair of 
simulations we can boost the statistical precision of the dark matter 
only simulations by factors of ~ 30 — 40, without biasing any of the 
results (Chuang et al. 2019). 

These simulation are complemented by full-physics hydrodynam- 
ical simulations with AREpo using the ILLustrisTNG galaxy forma- 
tion model (Weinberger et al. 2017; Pillepich et al. 2018a), employing 
two different box sizes: L = 500 h-'Mpc = 738.1 Mpc ~ 740 Mpc 
and L = 125h7!Mpe = 184.5Mpe ~ 185Mpc. The combined 
number of dark matter particles and gaseous cells for each box are 
2 x 43207 and 2 x 10807, respectively. These simulations match 
the mass resolution of the MTNG740 dark matter-only run, i.e., 
the corresponding dark matter and baryonic mass resolutions are 
1.12 x 108 h7!Mo and 2 x 107 h7!Mo. This lies intermediate to 
the mass resolution of the pairs TNG300/TNG300-2 and TNG100- 
2/TNG100-3, respectively, while the volume of MTNG740 is en- 
larged by a factor of 14.5 relative to TNG300. All the galaxy for- 
mation physics parameters of the MTNG runs are kept exactly the 
same as in the original ILLUsrrisTNG simulations, thereby allowing 
a direct assessment of numerical convergence by comparing galaxy 
properties to the TNG100 and even TNGSO series, and where ad- 
equate, allowing a Richardson extrapolating to a fiducial infinite 
resolution. There are however two significant changes to the physics 
model: Due to severe memory pressure in fitting the largest simu- 
lations onto the supercomputer available to us, we were forced to 
disable magnetic fields and to simplify the tracking of metallicity in 


the MTNG hydro simulations. This has, however, only a rather minor 
influence on galaxy properties. We refer to Pakmor et al. (2022) for 
a first in-depth analysis of the MTNG hydrodynamical simulations, 
and a comparison of its galaxy properties to TNG. 

We adopt the cosmological parameters given by Planck Collabora- 
tion (2016): Qm = Qegm+Qp = 0.3089, Q, = 0.0486, Qn = 0.6911, 
h = 0.6774, og = 0.8159 and ns = 0.9667 for this set of simula- 
tions, which have been used previously by ILLustrisTNG, such that 
the comparison to TNG is unaffected by any modification in cosmo- 
logical parameters. The initial conditions were generated at z = 63 
with second-order Lagrangian perturbation theory based on a new 
version of the NcENIC algorithm implemented into GADGET-4, us- 
ing the same linear theory power spectrum as used in the original 
ILLustRisTNG simulations. 

The high computational cost of the hydrodynamic simulation 
makes running a second realisation prohibitive, so we are here con- 
tent with constraining only the mode amplitudes, which already gives 
a good fraction of the benefit of the variance suppression technique, 
in particular for simple second-order statistics. The hydro simulation 
is matched to the “A” version of the corresponding dark matter only 
simulations. Note that the paired DM simulations in any case allow 
us to judge how important a second realisation with mirrored phases 
is for this technique. 

While our MTNG simulations constitute a transformative numeri- 
cal model for studying galaxy formation on the largest scales, it is not 
yet fully sufficient to address the cosmological science questions that 
require multi-Gpc? simulation volumes. To address this latter need, 
we boost the scope of our results with an extremely large N-body 
simulation, carried out with a two times better mass resolution than 
the original MILLENNIUM simulation (mp = 6.66 x 108 h-'!Mo) but 
with a much larger box size of 2.04 h-!Gpe = 3Gpe, plus the ad- 
dition of modelling of massive neutrinos, considering a sum of the 
neutrino masses equal to Xmy = 0.1eV. The dark matter particle 
number used in this simulation is 102403, augmented with a further 
2560° simulation particles used to represent the neutrinos. For the 
latter, we implemented a variant of the 6 method proposed by El- 
bers et al. (2021) in the Gapcet-4 code, which will be described in 
detail in Hernandez-Aguayo et al. (2022). 

The high computational cost of this run also precludes carrying 
out a further B version for the moment. However, we have done this 
for a set of three paired simulations of the same mass resolution but 
much smaller volume, and in which we have changed the neutrino 
mass, from Xm, = 0.3 eV, over Xm, = 0.1 eV, to Xm, = 0. For all 
of the neutrino simulations, we have also updated the cosmological 
model used to the newest DES-Y3 constraints (Abbott et al. 2022). 
Note that when we change the neutrino mass, we keep the total matter 
density at z = 0 constant. Also, we include photons and relativistic 
degrees of freedom in the background evolution, requiring changes 
in the GApGET-4 code, and a more sophisticated approach to set the 
initial conditions. Note, in particular, that in this case there is no 
closed form integration formula any more for computing the linear 
growth factor. 

The basic parameters and naming conventions of the MTNG sim- 
ulations are listed in Table 1. Analogous to ILLustrisTNG and for 
ease of comparison, we refer to our simulations with a name that 
is composed of the identifier “MTNG” followed directly by the box 
size in Mpc, slightly rounded where appropriate. Dark matter-only 
simulations additionally carry a designator “DM”, and simulations 
with neutrinos are labelled with the summed neutrino mass and a 
designator “v”. When different mass resolutions are available, they 
are distinguished at the end with a numerical identifier encoding the 
resolution level, with level “1” denoting the best available mass reso- 
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Type Run name Series Box size Neam Neas Ny Medm Meas/y x my €edm 
[h-'Mpc] [A'Mo] = [hk 'Mo] [eV] [A kpc] 

DM only MTNG740-DM(-1) A/B 500 43203 - - 1.32 x 108 - - 2.5 
MTNG740-DM-2 A/B 500 2160? - - 1.06 x 10° - - 5 
MTNG740-DM-3 A/B 500 10803 - - 8.50 x 10° - - 10 
MTNG740-DM-4 A/B 500 5404 - - 6.80 x 10!° - 7 20 
MTNG740-DM-5 A/B 500 270° - - 5.44x 10!! - - 40 
MTNG185-DM A 125 10803 - - 1.32 x 108 - - 2.5 

Hydro MTNG740 A 500 43207 ~—- 4320 - 1.12108 = 2.00x 107 - 2.5 
MTNG185 A 125 1080? ~—-1080 = 1.12x108  2.00x 107 = 2.5 
Neutrinos | MTNG3000-DM-0.1v A 2040 10240? = 2560 6.66x 108 =3.26 x 108 0.1 4 
MTNG1500-DM-0.1v A 1020 51203 = 12807 = 6.66 108 ~=— 3.26 x 108 0.1 4 
MTNG630-DM-0.3v A/B 430 2160% = 540° —s-6.54x 108 ~=— 9.76 x 108 0.3 4 
MTNG630-DM-0.1y A/B 430 21607 = 5407 6.66 x 108 ~— 3.26 x 108 0.1 4 
MTNG630-DM-0.0v A/B 430 2160° = 540°. «6.66 x 108 = 0.0 4 


Table 1. Specifications of the simulations of the MILLENNIUMTNG project introduced in this paper. The fiducial MTNG runs cover a volume of 5007 h~7Mpc? ~ 
740° Mpc? with resolution elements varying from 4320° (highest, level 1) to 270° (lowest, level 5), spaced by a factor of 8. Our naming convention uses the tag 
“MTNG” followed directly by the box size in Mpc, and optionally the identifier “DM” for dark matter-only runs followed by the resolution level, in analogy to the 
convention of ILLustrisTNG. Where needed for clarity, we append the letters “A” or “B” to distinguish the two different realisations run. The full hydrodynamical 
runs match the resolution of the MTNG740-DM simulations in two different volumes, 500° h~>Mpce? and 125° h~3Mpc?. In addition, we report the mass of 
dark matter particles and the initial mass of gaseous cells used in the full-physics runs, as well as the gravitational softening length. The last block of rows shows 
the specifications of the neutrino runs with box sizes of 2040 h~!Mpc, 1020 h~!Mpc and 430 h~!Mpe, with neutrinos of mass =m, = 100 meV (available for 
all box sizes), and &m,, = 300 meV and & m,, = 0 meV (for MTNG630 runs only). 


lution. When we need to explicitly distinguish the two different reali- 
sations, we append the letters “A” or “B” to the name. The simulations 
of the MILLENNIUMTNG project have been carried out for the most 
part on the SuperMUC-NG supercomputer at the Leibniz Computing 
Center. The large volume neutrino run (MTNG3000-DM-0.1yv) has 
been done on the Cosma8* machine at Durham University, while a 
few of the smaller simulations have been carried out on machines 
operated the Max Planck Computing and Data Facility (MPCDF)’. 


3 DATA PRODUCTS 


For the MILLENNIUMTNG project, we adopted a new outputting strat- 
egy designed to allow the construction of halo merger trees at high 
time resolution without the need to produce a large number of time- 
slices (traditionally called snapshots) on disk. This meant that much 
of the necessary halo finding and halo linking across time had to 
move from postprocessing to an on-the-fly treatment. Making this 
possible has been one of the main new features of the GapceET-4 
code, whose routines we use for this purpose. 

Furthermore, we wanted to simplify the comparison with deep 
wide-angle observational data by making use of a lightcone out- 
putting routine that detects crossings of simulation particles or cell 
trajectories with the past backwards lightcone of a fiducial observer 
position. We realise multiple lightcones of various depth and angular 
extent, and use them, in particular, to also create projected density 
shells for studies of weak gravitational lensing. Below, we give more 
details about the various simulation outputs produced by the MTNG 
simulations. 


2 nttps://www.dur.ac.uk/icc/cosma/cosma8 
3 https: //www.mpcdf.mpg.de/services/supercomput ing 
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3.1 Output times 


To produce determinations of group and subhalo catalogues, mea- 
surements of matter power spectra, and for (occasionally) storing 
snapshot particle information, the MTNG simulations define 265 
output times between redshift z = 30 and the present time, z = 0. 
The output spacing is constant in the logarithm of the scale factor, 
with three regions in which the spacing differs by a factor of 2 each, 
such that the resulting output frequency is finer at low redshift than 
at very high redshift, as follows: 


e Alog(a) = 0.0325 for 
e Alog(a) = 0.0162 for 
e Alog(a) = 0.0081 for 


10 < z< 30 (32 times). 
3<z<10 (62 times). 
O<z<3 (171 times). 


For the neutrino runs, we have reduced the output frequency by 
a factor of two to arrive at just 133 output times, omitting every 
second of the output times defined above. This was done in order 
to save some computation time, and especially disk storage space, 
and has also been motivated by our finding that the resulting time 
resolution for the updating of the group catalogue is sufficient for the 
semi-analytic galaxy formation code. 


3.2 Group and subhalo catalogues 


At all output times defined above, our simulation codes GADGET- 
4 and Arepo, respectively, have run the Friends-of-Friends (FoF) 
group finding algorithm (Davis et al. 1985) combined with the 
SuBFIND-HBT substructure finder to compute group and subhalo 
catalogues on-the-fly. The corresponding algorithms are described 
in the Gapcet-4 paper (Springel et al. 2021). Compared to the tra- 
ditional SuBrinp algorithm, the ‘hierachical bound tracing’ (Han 
et al. 2018, HBT) extension yields better tracking of subhaloes espe- 
cially close to pericentre. For each FoF group, we compute spherical 
overdensity mass estimates around the particle with the minimum 
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Figure 1. Projected dark matter density in one of the lightcones of the MTNG740 dark matter-only simulations (here the MTNG740-DM-2-A simulation was 
selected, for definiteness). The observer at the present time (z = 0) is located at the bottom centre of the figure, while the lightcone extents out to z = 2. The dark 
matter particles are displayed using comoving coordinates, and the comoving thickness of the projected inclined slice is 15 h~!Mpc = 22.14 Mpc, independent 
of distance. We also show two nested zoom-in regions (spherical insets) that illustrate the well-known spider web-like large-scale structure and a dense region 
around a forming galaxy cluster. These zoomed inset images have diameters of 400 Mpc and 40 Mpc, respectively. 


gravitational potential* which is taken as centre of the group. For 
each gravitationally bound subhalo, a number of further properties 
are computed, among them the maximum circular velocity, a mea- 
sure of the environmental density, shape information, and further 
properties. 

The particle [Ds that make up each group or subhalo are not 
explicitly stored. Instead, the simulation codes remember this mem- 
bership information until the next group catalogue is computed, and 
then use it to determine descendant and progenitor pointers that link 
two subsequent group catalogues. These pointers are stored on disk 
alongside the group and subhalo catalogues themselves, and consti- 
tute the basic information needed to construct detailed merger trees 
in a postprocessing step detailed below. 


3.3 Snapshots 


The above procedure obviates the need to store full snapshot data 
for all output times desired in a high time-resolution merger tree. 
In fact, none are in principle needed. However, for certain analysis 
(including unforeseen ones) full time slices are still needed, but then 
only a comparatively small number is usually sufficient, making still 
an order of magnitude reduction of the data volume possible. We 
have decided to store only 10 full particle snapshots at redshifts 


4 Computed for the set of particles making up the group. 


z= 0, 0.25, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, and 7.0>. For the neutrino 
simulations, we reduced this even further to just 5 snapshot times. 

For output times where a corresponding full snapshot is saved, we 
note that the particle data is stored in an ordered and nested fashion, 
with the largest FoF halo coming first (with the others following 
in descending order), with the particles making up each subhalo 
within a FoF group being stored in descending order of subhalo size 
as well. Finally, within each subhalo, the particles are ordered by 
increasing binding energy, so that the most-bound particle comes 
first in each substructure. This storage scheme thus makes it possible 
to selectively load the particles making up each group or subhalo, 
as they are stored consecutively on disk with a starting offset that is 
known beforehand. Also, this does not require a separate storage of 
the IDs that make up a certain group or subhalo. 

For the special application of semi-analytic galaxy formation, we 
have actually produced an additional set of snapshot files at all output 
times. These contain only those particles that have been a most- 
bound particle of a subhalo sometime in the past. These particles 
can be used to approximately track in semi-analytic models so-called 
orphaned galaxies whose dark matter substructures are disrupted by 
tidal forces before the corresponding galaxies are predicted to have 
merged with their central galaxies. The fraction of these particles 
grows monotonically in time and reaches a few percent by z = 0. This 
means that even with 265 output times as used here, the cumulative 


5 The corresponding snapshot numbers for the main MTNG runs are 264, 
237, 214, 179, 151, 129, 094, 080, 069, and 051. 
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Figure 2. Particle mass resolution impact on the non-linear power spectrum of the MTNG740-DM runs at z = 0 (left panel) and z = 1 (right panel). In each 
case we show the average of the A and B realisations that we carried out. Different colours represent the different number of particles, as labelled, corresponding 
to our resolution levels 1 — 5. In both panels, the vertical dashed line indicates the fundamental mode of the box, and the diagonal solid grey lines give the 
shot-noise contribution for each case. The solid black lines display the linear theory prediction. The lower subpanels show the relative difference with respect to 
the highest resolution MTNG740-DM runs, and the horizontal dashed lines mark the 1% difference interval. The dashed lines in the lower subpanels correspond 
to the relative difference of the PS measurements without subtracting the shot-noise contribution. 


data volume of these special most-bound particle snapshots is still 
small. 


3.4 Power spectra measurements 


Matter power spectra are measured for all defined output times (both 
for the total matter, and, where available, separately for baryons, 
dark matter, and neutrinos). Three power spectra measurements are 
done at each output time, applying the folding technique (Jenkins 
et al. 1998) twice with folding factors of 16 and 162. This allows 
the measurements to be combined such that they yield a coverage 
of the full spatial dynamic range of the simulations, up to k—ranges 
that reach well below the gravitational softening scale. The power 
spectra are output in a finely binned fashion such that they can be 
conveniently band-averaged for a new binning, if desired. 


3.5 Particle lightcones 


For each simulation, we produce a set of different lightcones for a 
fiducial observer located at the origin of the computational box®. 
The periodic simulation box is replicated automatically to the extent 
necessary to fill the geometry of the specified lightcone. We have 
created five different lightcones with the following definitions: 


e Cone 0: full-sky particle lightcone between z = 0—0.4, extend- 
ing to a comoving distance ~ 1090 h~!Mpe. 

e Cone 1: particle lightcone covering one octant on the sky be- 
tween z = 0 — 1.5, reaching comoving distance ~ 3050 h-'Mpe. 


® Note that due to the translational symmetry induced by our periodic bound- 
ary conditions, this point is not special in any way. 
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e Cone 2: a pencil beam particle lightcone with a square-shaped 
footprint of area (10 deg)* at an oblique angle, between z = 0 — 5, 
hence extending out to a comoving distance ~ 5390 A~!Mpc. 

e Cone 3: a disk-like particle lightcone with a comoving thickness 
of 15 h-'Mpc (tilted against the principal coordinate planes), over 
the redshift z = 0 — 2, and thus extending to a comoving distance 
~ 3600 h7!Mpe. 

e Cone 4: a full sky lightcone between z = 0 — 5, but only con- 
taining ‘most-bound’ particles as described for the partial snapshot, 
yielding a comoving distance ~ 5390 h-'Mpe. 


Note that while Cone | and 2 redundantly comprise some of the 
particle data that is contained in Cone 0, they go out to considerably 
deeper redshift, which is the reason they were added in the first place. 
The primary purpose of Cone 3 is to allow a two-dimensional visu- 
alisation and analysis of structure out to higher redshift than possible 
by extracting this data from a full-sky lightcone. Finally, Cone 4 can 
be used in semi-analytic modelling of galaxy formation to improve 
the orbital treatment of galaxies by allowing a precise determination 
of when subhaloes or orphans (whose position is both marked by 
particles which have been, or still are, a most-bound particle of a 
subhalo) cross the past-backwards lightcone, thus yielding accurate 
phase-space information for these events. We note in passing that for 
the neutrino runs, we added a second pencil beam lightcone pointing 
into a different direction than Cone 2. 

As an illustration of the lightcone outputs, in Fig. 1 we show a cut- 
out from Cone 3 (the disk-like lightcone). The observer is located at 
z = 0 at the bottom centre of the figure, while the lightcone extends 
out to z = 2, which is reached at the outer perimeter of the displayed 
half-sphere. The right circular inset shows a zoomed region of the 
lightcone, centered on a forming galaxy cluster. This highlights the 
filaments, knots and voids that comprise the large-scale structure of 
the Universe. A further zoom by another factor of 10 (left spherical 
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Figure 3. Power spectrum measurements for our MTNG3000-DM-0.1v run at the initial time (left panel, z = 63) and the final time (right panel, z = 0). In both 
cases we show the measured power spectrum of the cold matter particles (CDM and baryons; solid blue line) and of the neutrino particles (solid red line), and 
we compare to the expected linear theory power spectrum (thin dotted curves) as computed by CAMB (Lewis et al. 2000) for this cosmology (which includes 
neutrinos that transition from the relativistic into the non-relativistic regime, as well as a photon background). The dashed vertical line shows the fundamental 
mode of the box, while the vertical thin dotted line gives the Nyquist frequency of the initial DM particle grid. The dotted horizontal lines gives the nominal 
shot noise of the two particle sets, 10240? for the DM, and 2560° for the neutrinos. Despite this, the actual shot noise realised by the neutrinos at late times 
is much lower, thanks to the 6 f -simulation technique (Elbers et al. 2021), and is in fact close to the one of the dark matter at late times. In the bottom panel, 
we give the relative deviation of the measured power spectrum modes relative to linear theory. Note that in the initial conditions, we deliberately create a small 
boost in the linear power at the largest scales in order to offset the fact that our code treats super-horizon modes with Newtonian gravity. 


inset with radius 20 Mpc) displays the distribution of individual dark 
matter haloes and their embedded substructures in the dense region 
of the protocluster. 


3.6 Mass-shell outputs 


For weak gravitational lensing applications, we additionally create 
onion-like shells with projections of a fiducial full-sky particle light- 
cone (which itself is not output to disk due to its prohibitive size) 
onto a healpix tesselation of the sky. The comoving depth of these 
shells is 25 h~'Mpe, and they go out to redshift z = 5, giving 216 
such maps in total. The number of equal-area pixels in each map 
is Npix = 12 Ne der Where Ngide is the resolution parameter of the 
healpix tesselation algorithm. For our highest resolution maps we go 
up to Neide = 12288, yielding 1.8 billion pixels and a 0.28 arcmin 
angular resolution of these mass maps. Each of the pixels simply 
contains the total mass of all the particles that fall within the corre- 
sponding solid area of the pixel. 


3.7 Merger trees 


We build merger trees for all the simulations primarily in order to con- 
struct mock galaxy catalogues using semi-analytical models. In prin- 
ciple, the group/subhalo catalogues and the descendant/progenitor 
information produced on-the-fly during the simulation runs already 
contain all the data needed for the merger tree. However, to efficiently 
work with this data (in particular to avoid excessive I/O times to col- 
lect it from many different files), it is prudent to rearrange the data 
such that the subhaloes that are linked together in a single tree are 


also stored together. This final step in the merger tree construction is 
carried out in a postprocessing step using the methods implemented 
in the Gapcet-4 code (Springel et al. 2021). 

The result are tree catalogues where each tree is self-contained 
in the sense that all progenitor and descendant pointers only lead to 
other subhaloes contained in the same tree. Also, all subhaloes in a 
common FoF group are always in the same tree. A single tree is thus 
sufficient for running semi-analytic models of galaxy formation such 
as L-Gavaxies (Henriques et al. 2015), which can therefore process 
the trees of a simulation in an “embarrasingly parallel” fashion. As 
a further convenient data product, GApGEt-4 creates auxiliary files 
that tell for each subhalo in which tree, and at which place within 
the tree, the corresponding subhalo can be found. One can then 
selectively load this tree, if desired, to examine, for example, the 
history and fate of the chosen subhalo. 


4 PRECISION PREDICTIONS FOR MATTER AND HALO 
STATISTICS 


The MILLENNIUMTNG simulations aim to assist ongoing and future 
galaxy surveys (e.g., DESI, Euclid and PFS) by making accurate 
predictions of the clustering of matter over cosmic time (and thus the 
clustering of different types of galaxies). In this section we verify the 
statistical power of the MTNG simulations by performing a series of 
convergence tests, by highlighting the advantage of running fixed- 
and-paired simulations, and by assessing the accuracy of the matter 
and halo statistics predicted by published emulators. In the following, 
we will always show the mean of the measurements for the A and B 
realisations where available, unless stated otherwise. 
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4.1 Numerical convergence assessment of matter and halo 
statistics 


4.1.1 Non-linear matter power spectra 


The non-linear matter power spectrum is a key outcome of struc- 
ture formation and allows us to test the convergence of numerical 
simulations, and to establish the smallest scales at which accurate re- 
sults can be obtained (see e.g. Heitmann et al. 2010; Schneider et al. 
2016; Grove et al. 2022). Furthermore, large galaxy surveys require 
accurate predictions of the matter power spectrum in the nonlinear 
regime, i.e. up to scales k ~ 10h Mpc7! with an accuracy of order 
1 per cent (LSST Science Collaboration 2009; Laureijs et al. 2011; 
DESI Collaboration 2016). 

As described earlier, we have measured the matter power spectra 
for our MTNG simulations for all defined output times (both for the 
total matter, and where available, separately for baryons, dark matter, 
and neutrinos), and with the folding technique down to k—ranges that 
reach well below the softening scale, thus covering the full spatial 
dynamic range of the simulations. 

In order to illustrate the convergence of our MTNG simulations, 
Figure 2 displays the measured (dimensionless) non-linear matter 
power spectrum from our MTNG740-DM runs at z = 0 (left panel) 
and z = 1 (right panel). Each coloured curve shows the average of 
the A and B realisations for each resolution, while the black solid 
line displays the linear theory prediction. The dimensionless power 
spectrum is expressed as, 


2 k3 
A*(k) = aya mtk) » (1) 


where P,,(k) is the non-linear matter power spectrum measured 
from our MTNG simulations. The vertical dashed line represents the 
fundamental mode of the box, given by 


2n 
Kpox = > 2 
box L (2) 


where L = 500h~!Mpc is the box length of the MTNG740-DM 
simulations. In addition, the diagonal lines display the Poisson shot- 
noise contribution, 


k3 
A? (k) shot = Fa P shot , (3) 


with Pgnot = L3/ Np and Np being the total number of particles. We 
have here subtracted the shot-noise from our power spectra measure- 
ments (but see Maleubre et al. 2022). 

From the upper panels of Fig. 2 we can clearly see that the size 
of the MTNG740 simulations is enough to measure the clustering 
of matter up to the BAO scale, as explored further below. Also, we 
have a first impression of the mass resolution impact on the power 
spectra measurements, low-resolution simulations predict a lack of 
power on small-scales compared with higher resolution runs. We 
quantify the dependence of mass resolution on the power spectrum 
through the relative difference between the lower resolution runs of 
level-5 to 2 (denoted based on their particle numbers as N270, N540, 
N1080 and N2160) and the highest resolution level-1 simulation 
(N4320, see the lower subpanels of Fig. 2). We find sub percent 
(< 1%) agreement at z =Oonscalesk <0.2h Mpc7!, 0.6h Mpce™!, 
2hMpc™! and 104 Mpc™!, for the N270, N540, N1080 and N2160 
simulations, respectively. In addition, the dashed lines in the lower 
subpanels of Fig. 2 indicate the relative difference between the PS 
measurements without subtracting the shot-noise contribution. We 
find a slightly improved agreement between the lower resolution 
and the level-1 runs on small-scales. Similar values are obtained at 
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z = 1, as seen in the right panel of Fig. 2. Based on these results, we 
conclude that cosmological simulations with mass resolution better 
than mp = 1.06 x 10° h-'!Mo are needed to achieve the precision of 
theoretical predictions required by the current and upcoming galaxy 
surveys. 


To illustrate the accuracy of our MTNG3000-DM-0. 1 v simulation, 
we show in Figure 3 the power spectrum measurements for the cold 
(dark) matter (CDM and baryons) and neutrino components both at 
the initial time, z = 63, and at the final time z = 0. We compare 
directly to the linear theory prediction of CAMB (Lewis et al. 2000) 
for the corresponding cosmology, assuming that baryons are present 
but represented by the dark matter as well. The neutrinos consist of 
two degenerate 50 meV species, and a third massless neutrino. We 
also include the photon background, and have modified GADGET- 
4 such that it correctly accounts for relativistic backgrounds in the 
clustering evolution. At the initial time, we however deliberately put 
in a small correction in the form of a power increase of up to ~ 4% 
at the fundamental mode, because at this time these modes are still 
outside of the horizon, meaning that our Newtonian code would not 
capture their faster growth accurately. The correction is computed 
such that at low redshift our results match linear theory precisely. 
The bottom panel of the z = 0 result shows that this succeeds with 
impressive precision. Our simulation reproduces the expected linear 
theory results on the largest scales to fractions of a percent, both in 
the dark matter and in the neutrino components. 


Note that the power in the cold dark matter necessarily cuts off 
at the Nyquist frequency in the initial conditions, at which point 
the cold dark matter power actually lies far below the shot noise 
limit. The ‘coldness’ of the dark matter allows this sub-Poissonian 
distribution to be preserved during the subsequent evolution. As our 
power spectrum measurement extends to much smaller scales, we 
find, however, for scales below the Nyquist frequency in the initial 
conditions the shot noise power, simply because these scales are 
undetermined in the initial conditions. They are however filled in 
during non-linear evolution by power transfer from larger scales, 
such that the non-linear power at z = 0 can be predicted accurately 
well below the initial Nyquist frequency. Eventually, however, the 
dark matter power spectrum hits the shot noise limit in the deeply 
non-linear regime. 


For the neutrinos, the initially imprinted power spectrum would 
normally not be preserved below the neutrino shot noise, due to their 
extreme ‘hotness’, letting them move with velocities close to the speed 
of light. Interestingly, this unfavourable outcome does not occur for 
the neutrinos, and their effective shot noise (horizontal plateau of 
the measurement) stays at a level far below the nominally expected 
shot-noise for the 2560° neutrino particles we used. This is due to the 
6 f-method (Elbers et al. 2021) we employed to simulate the neutrino 
component, which allows us to keep the shot-noise level at late times 
nearly a factor of 100 lower than naively expected, making it in fact 
very close to that of the dark matter particles, a situation we intended 
to achieve with our choice of the particle numbers. 


In the remainder of this paper, we shall focus on an analysis of the 
500 h-'Mpc simulations (MTNG740), deferring a further analysis 
of the neutrino predictions to forthcoming work (Hernandez-Aguayo 
et al. 2022). We proceed in this way because we consider it prudent 
to first validate the foundations of the MTNG simulation project 
and test for numerical convergence and the impact of baryons on 
basic clustering statistics, before exploring the additional effects of 
neutrinos in more quantitive detail. 
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Figure 4. Convergence tests of the measured halo mass (M09,.) function from the MTNG740-DM runs at z = 0 (left panel) and z = 1 (right panel). The lower 
subpanels show the relative difference with respect to the highest resolution MTNG740-DM simulations, with the horizontal dashed lines representing a 1% 
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Figure 5. Similar to Fig. 4, but for the measured halo mass function using the Friends-of-Friends halo mass (Mfop) definition. The solid and dashed lines 
correspond to the measured HMF with and without the Mfor correction given by Eq. (4), respectively. 


4.1.2 Halo mass function and halo clustering 


Dark matter haloes are the building blocks of the cosmic structures 
since they host the different types of galaxies in the Universe. There- 
fore, the study of their abundance and clustering is important for 
modern cosmological analysis. The abundance of dark matter haloes 
can be quantified by the halo mass function (HMF), which gives the 
number of dark matter haloes as a function of halo mass, cosmic 
time and cosmological parameters, in different mass intervals. Since 
cosmic structure formation is a hierarchical process, small haloes 
form first and their subsequent mergers build ever larger groups, up 


to haloes with the mass of galaxy clusters, which are the most mas- 
sive objects in the Universe. We expect to find a higher abundance of 
less massive objects at early-times (or higher redshifts), whereas the 
most massive haloes are rare (and hence require large simulation vol- 
umes for good statistics) and form comparatively late. In any case, 
the HMF can be considered one of the most important quantities 
predicted by cosmological simulations. Exploring the completeness 
and convergence of the HMF, which is strongly related to the mass 
resolution and the simulated volume, is thus of great importance to 
accurately predict galaxy statistics. 
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Figure 6. Convergence test of the real-space halo two-point correlation functions at z = 0 (left panel) and z = 1 (right panel), for three different halo samples 
with masses Mog9- > 10!!-5 h-!Mo (solid lines), Moog > 10!25 h7!Mo (dashed lines) and Mzo9¢ > 10!3-> h-!Mo (dash-dotted lines), based on the 
MTNG740-DM simulations. The relative difference between the low and the highest resolution simulations is shown in the lower subpanels. The different 
resolution levels are distinguished by their particle number per dimension, as labelled. 
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Figure 7. Box size impact on the measured non-linear matter power spectrum (left panel) and the halo mass function (right panel) at z = 0 (blue lines), z = 1 
(orange lines), z = 2 (green lines), and z = 3 (red lines) from the MTNG740-DM-1 (dashed lines) and MTNG185-DM (solid lines) simulations, the dark 
matter-only versions of the MTNG hydro simulations. The lower subpanels show the relative difference with respect to the larger MTNG740-DM- 1 simulation, 
where the horizontal dashed lines represent differences of 5% (left panel) and 1% (right panel), respectively. Note that for MTNG740-DM we use the average 
result of two paired realizations, whereas for MTNG185-DM only one (but still ‘fixed’) realization is available. 


In Figure 4, we show the measured halo mass functions at z = 0 
(left panel) and z = 1 (right panel) from the MTNG740-DM dark 
matter-only simulations. We use here the Mzgo- spherical overden- 
sity halo mass definition which corresponds to the mass enclosed 
within a radius in which the density is 200 times the critical density 
of the Universe. From the upper panels of Fig. 4, we see that for differ- 
ent resolutions we can find very similar HMFs at the high-mass end. 
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In addition, we note a cutoff at the low-mass end due to the limited 
number of particles in low-mass haloes. This cutoff moves system- 
atically to high-mass haloes when we reduce the resolution of the 
simulations. The cutoff scale can be used to assess the completness 
of the HMFs; i.e., the point where the HMFs show the largest differ- 
ences due to the resolution of the simulations. The level-2 (N2160) 
and level-1 (N4320) simulations predict almost the same number of 
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haloes for masses Moog¢ > 10!2 h7'Mo at both z = Oand z = 1. Sub 
percent agreement between the N1080 and N4320 runs is found for 
haloes more massive than M9, = 10!3 A~'!Mo. The very low res- 
olution of the N540 and N270 runs have, however, a greater impact 
on their measured HMFs, as we can see through the large differences 
when compared with the HMF from the N4320 run. Furthermore, 
it appears that the HMFs for the N540 and N270 simulations is not 
converged at z = 1 for any halo mass, since they fall short of produc- 
ing the same number of massive haloes as the N4320 run (see right 
panel of Fig. 4). Nevertheless, these low-resolution simulations are 
still useful to test and debug analysis pipelines. 

As an additional convergence test of the halo population, Fig. 5 dis- 
plays the measured HMFs based on the Friends-of-Friends (Mpfor) 
halo mass definition at redshifts z = 0 and z = 1. However, it has been 
shown that the FoF definition has discreteness limitations when esti- 
mating the halo mass from low-resolution simulations, producing on 
average an overestimation of the HMF measured from low-resolution 
simulations (Warren et al. 2006). We observe this trend when mea- 
suring the HMF directly from our halo (FoF group) catalogues as 
shown in the lower subpanels of Fig. 5 (dashed lines). We find an 
excess of haloes between 2 — 10 percent when comparing the level 
2 — 5 runs with our flagship MTNG740-DM-1 simulations. For this 
reason we apply a correction suggested by Warren et al. (2006) to 
each FoF halo mass, given by 


Mror = (1—- N50) MS, (4) 


where Np is the number of member particles of a halo and maim is 
the FoF mass of each halo in the original catalogue. 

We then find a ~ 1 per cent agreement between the FoF halo 
mass functions of the N2160 and N4320 runs for haloes with masses 
Mror > 10!!-5 h-!Mo, at both redshifts. We still see a slight en- 
hancement of the FoF mass for the N1080, N540 and N270 runs 
compared to the higher resolution simulations, amounting to 3 — 5 
per cent at the present time, but these differences are smaller at early 
times, e.g., at z = 1, where the differences are 2% (see the right 
lower-subpanel of Fig. 5). 

We next explore the real-space halo two-point correlation function 
and consider its numerical convergence for our simulations. The 
correlation function gives the excess probability to find a pair of 
tracers (e.g. dark matter haloes, or galaxies) separated by a distance r 
compared to a random, uniform distribution. Since haloes are biased 
tracers of the underlying dark matter field, their clustering amplitude 
can be different for different halo populations, e.g., massive haloes 
are more strongly clustered than low-mass haloes (the linear halo 
bias will be discussed in Sec. 5.2). 

We therefore measure the real-space clustering for three dif- 
ferent halo samples with fixed Mzo9- halo mass, and using 20 
logarithmically-spaced radial bins between 0.5 < r/[ h~'!Mpc] < 50 
using the NBopyxir toolkit (Hand et al. 2018). We select haloes with 
masses Mo99¢ > 10!!° A-!Mo (sample 1), Maoo¢ > 10!23 A-'Mo 
(sample 2), and Mrq9¢ > 10!3-> h-!Mo (sample 3). Note that not all 
simulations are able to resolve haloes for each halo sample (see the 
completeness of the Mog9- HMF in Fig. 4), for this reason sample 1 
contains haloes from the N4320 and N2160 runs only, sample 2 cov- 
ers the N4320, N2160 and N1080 runs, and haloes of sample 3 can 
be found in all of the MTNG runs (N4320, N2160, N1080, N540 and 
N270). 

The measured halo clustering is shown in the upper subpanels 
of Fig. 6 at redshifts z = 0 and 1. We can clearly see that the most 
massive haloes (sample 3; dash-dotted lines) have a higher clustering 
amplitude, which means that they are more biased than low-mass 
haloes (sample | and 2; solid and dashed lines, respectively). In the 


lower subpanels of Fig. 6, we compare the halo clustering of the 
different resolution runs with respect to N4320. The halo correlation 
functions of the N4320 and N2160 runs agree within | per cent for all 
the samples at z = 0 and at all scales, however, the same agreement is 
found at z = | for all length-scales r > 1 h~!Mpe. The clustering of 
samples 2 and 3 of the N1080 case show a slightly larger difference 
of ~ 3% in comparison with the N4320 runs for scales greater than 
1 h~'Mpc; similar differences are found for sample 3 of the N540 
simulation (see orange dash-dotted lines). The lowest resolution runs 
(N270) produce a larger difference of ~ 6 per cent, this is due to the 
incompletness of the halo mass function displayed in Fig. 4, where 
we found a deficit of more than 10 per cent for haloes with masses 
Mo90¢ > 10!3° h7 Mo. 


4.1.3 Box size effects 


Finally, we show the impact of the box-size on the measured matter 
power spectrum and the halo mass function of the MTNG simulations 
in the left and right panel of Figure 7, respectively. It is important 
to perform convergence tests by varying the size of the cosmologi- 
cal box and keeping the same particle resolution, to understand the 
limitations in cosmological calculations due to the finite volume of 
numerical simulations. To this end we compare the measurements 
for the MTNG740-DM and MTNG185-DM simulations, which are 
the dark matter-only counterparts of the corresponding full-physics 
simulations. The MTNG185 run has a volume 64 times smaller than 
our flagship MTNG740 run, while keeping the same mass resolution 
(see Table 1). 

The box-size effects on the non-linear power spectrum at z = 0, 1, 
2 and 3 are displayed in the left panel of Fig. 7. The vertical dashed 
lines show the fundamental mode of each box given by Eq. (2), 
which delimits the minimum scale we can measure the power spectra 
from the corresponding simulations, i.e., kt,599 ~ 0.0125 h Mpc7! 
and ky125 ~ 0.054 Mpc™!, respectively. The clustering of matter 
agrees very well between the two boxes at all redshifts; we find only 
a ~ 2% difference on scales k = 0.13h Mpc! in all cases (see the 
left lower-subpanel of Fig. 7). We observe larger deviations of up to 
20% for the large-scale modes close to the k,125 value. This is due 
to the limited box size of the MTNG185 simulation, which does not 
allow us to measure the clustering on these large scales with sufficient 
precision, both because they are already affected by mild non-linear 
evolution (which is distorted by missing, still larger scale modes) and 
because we have run only one MTNG185-DM realization, which is 
thus still dominated more strongly by sample variance as we can not 
eliminate in this case the leading high-order perturbations around 
linear evolution by means of the pairing technique. 

The evolution of the HMF from z = 3 to z = 0 of the MTNG740- 
DM and MTNG185-DM boxes is shown in the right panel of Fig. 7. 
At first glance, we can see that the HMFs of the MTNG185 run 
follow the shapes of their MTNG740 counterparts very well, showing 
completness down to the same value of My99- = 10!9 h-!Mo, which 
corresponds to well-resolved dark matter haloes with 75 particles. 
The agreement between the HMFs is much better than 1% for haloes 
with masses Myg9- < 10!3 h7!Mo at all redshifts. We find some 
small fluctuations at the high-mass end (Mago. > 10/3 h-'Mo) 
due to the limited volume of the MTNG185 simulation, which does 
not contain enough massive haloes, such as groups and clusters-like 
objects, to avoid being severely affected by small number statistics. 

Thus far, we have demonstrated that our MTNG740-DM runs 
reach enough cosmological volume and mass resolution to show 
good convergence in the non-linear matter power spectrum and halo 
mass functions. This should make especially the level-1 (N4320) 
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Figure 8. The ratio between the measured non-linear matter power spectrum of the MTNG740-DM-1-A (red lines) and -B (blue lines) simulations relative to 
the linear theory power spectrum prediction, at different redshifts, as labelled. The green line represents the mean of the A and B realisations. 
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Figure 9. Comparison of the measured power spectrum (left panel) and halo mass function (right panel) between the mean over 100 independent Gaussian 
realisations (dashed lines) and the A— and B- fixed-and-paired (solid lines) MTNG740-DM-5 runs at z = 0 (blue lines), z = 0.5 (green lines), z = | (red lines), 
and z = 2 (cyan lines). The lower subpanels show the relative difference between the results from the 100 Gaussian realisation (R-mean)s and the A/B pair, 
while the dashed lines indicate a 1% (left panel) and 2% (right panel) difference. 


and level-2 (N2160) runs very helpful for the generation of mock 
catalogues for large-scale galaxy surveys such as DESI, Euclid or 
PFS. 
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4.2 Reducing cosmic variance with fixed-and-paired 
realisations 


As mentioned above, we have used the fixed-and-paired technique 
by Angulo & Pontzen (2016) to produce two DM-only cosmological 
simulations for a given model. The method consists of running two 
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realisations of the initial density perturbation fields, both set-up with 
Fourier modes amplitudes fixed to the ensemble average power spec- 
trum, and in addition, both are paired through the use of opposite 
phases for each mode (i.e., with a phase difference of 2, or equiva- 
lently, with a sign reversal of the density perturbation field). For pure 
linear theory, the pegging of the mode amplitudes to x »/P(k) repro- 
duces the linear theory input power spectrum by construction at all 
times, even when considering each of these realisations individually. 
However, mild non-linear evolution will still introduce deviations of 
the mean power in the large-scale modes. But, as Angulo & Pontzen 
(2016) show, these can be cancelled to leading order when the re- 
sults of the two realisations are averaged. In this way, such a pair of 
simulations can yield equivalent results as obtained from the mean 
of many independent realisations. 

Figure 8 shows the evolution of the ratio between the measured 
non-linear dark matter power spectrum from the MTNG740-DM- 1- 
A and -B realisations and the linear theory prediction, from z = 30 to 
z = 0. We can see that taking the average of the two fixed-and-paired 
realisations (green lines) eliminates the higher order fluctuations (red 
and blue lines) around linear theory, allowing us to obtain much more 
accurate cosmological predictions at large scales for these still com- 
paratively ‘limited-volume’ simulations. Importantly, this method 
works well down to the present time on scales not yet strongly af- 
fected by non-linear evolution, even though individual modes already 
differ on these scales by 10s of percent from pure linear theory due 
to higher order terms affecting the evolution. 

As an additional consistency check of the advantage of running 
fixed-and-paired simulations to suppress the cosmic variance, we 
have run 100 independent, fully Gaussian realisations with the same 
specifications as the MTNG740-DM-5 simulations (i.e. box size, 
number of particles, and cosmological parameters kept the same). 
We then compare the matter clustering and halo mass functions 
of our fixed-and-paired simulations with those from the Gaussian 
realisations. 

The matter clustering results at z = 0, 0.5, 1 and 2 are shown in 
the left panel of Fig. 9. From the lower subpanel we can see that the 
average of the A and B realizations of MTNG740-DM-5, and the 
mean over the 100 independent realisations agree to 1% at scales 
k < kny, where kny is the Nyquist frequency (indicated by the 
vertical dashed line) given by, 


aN grid 
L 


ky = , (5) 


where Ngyiq is the number of particles per dimension, in this case 
Nerid = 270. Even on scales smaller than the Nyquist frequency, 
where power is created due to non-linear clustering, there is no 
systematic difference, only the random fluctuations become larger 
and reach levels of several percent. 

The halo population is also affected by cosmic variance, since 
haloes form and evolve from peaks in the dark matter density fluc- 
tuations. In order to quantify the impact of the variance suppression 
technique on this statistic, we investigate the differential halo mass 
functions at the same output times used above, at z = 0, 0.5, 1 and 
2, and give the result in the right panel of Fig. 9. The HMF shows 
qualitatively the same level of agreement as the power spectrum. The 
predicted halo population of one set of paired-and-fixed simulations 
is almost identical as the one obtained by averaging 100 independent 
realisations. The somewhat larger differences at the high-mass end 
are here due to the poor counting statistics in the exponential tail of 
the mass function that is invariably still present as a result of the still 
limited volume of the MTNG740 runs. 


4.3 Comparison with emulators 


It is well known that cosmological simulations are computationally 
quite expensive, especially if one wants to push towards the large 
volumes and the high-resolution required by modern galaxy surveys, 
as we have done here with MTNG. In order to alleviate the com- 
putational cost of running such numerical calculations, researchers 
have developed cosmological emulators based on fits to simulation 
results, with the goal to use them as interpolation (or even extrapo- 
lation) tools to predict the outcome of nonlinear structure formation 
for a variety of cosmological models in a very rapid, but ultimately 
still approximative fashion (see e.g., Heitmann et al. 2006; Habib 
et al. 2007). 

Here, we assess the accuracy of cosmological emulators for the 
non-linear matter power spectrum and the halo mass function when 
compared with our MTNG simulations. On the one hand, we use pre- 
dictions of the matter power spectrum from the CosmicEmu (Moran 
et al. 2022), Bacco (Angulo et al. 2021) and EucLipEMULATOR2 
(Knabenhans et al. 2021) emulators. All of these emulators were cal- 
ibrated using a large number of N-body simulations (O(100)) with 
different cosmological parameters’. 

Instead of comparing the full-shape of the power spectrum, we use 
the non-linear boost, which is the ratio between the non-linear matter 
clustering and its linear theory prediction, 


= Pmlk) 
cae Piin(k) 


The comparison of the boost factor, B(k), at z = Oand z = 1, isshown 
in the left and right panel of Fig. 10, respectively. At the present 
time (z = 0; lower-left panel), the predictions of the MTNG740-DM 
simulations and the emulators are consistent within 1% accuracy 
fork > 0.2h Mpc™!, with some fluctuations seen at the large-scale 
modes, but the differences are well within the 5% limit. At early 
times (z = 1; lower-right panel), we find slightly larger differences 
than 1% at the high-k modes, but this is consistent with the claimed 
accuracy of the emulators. 

Making fast predictions for the halo mass function is an important 
goal for emulators as well. In the left panel of Fig. 11, we compare 
the differential halo mass function by Bocquet et al. (2020) (which is 
a part of the CosmicEmu project) and our MTNG results at z = 
0 (blue lines) and z = 1 (orange lines). The HMF CosmicEMu 
was built using the Mago. halo mass definition and is limited to 
the mass range 10!3 < Myg9¢/[h~'Mo] < 10!9° at z = 0. The 
CosmicEmu haloes were identified with the FoF algorithm with a 
linking length parameter of b = 0.168, and a spherical overdensity 
(Mo09-) halo catalogue was built using the potential minimum of 
each FoF group as halo centre, in a similar way to the SUBFIND 
code. We find an agreement better than 5% for haloes with mass 
Mo9¢ < 10!5 A-!Mo (Map9¢ < 10!4h7!Mo) at z = 0 (z = 1). 
There are larger differences of ~ 10% at z = 1 at the high-mass end. 

We also make use of the DARKQuEstT emulator (Nishimichi et al. 
2019) to test the accuracy of the HMF. This HMF emulator was 
calibrated using the Myo, halo mass definition, which corresponds 
to the mass within a sphere with overdensity 200 times the mean 
background density of the Universe. Note that DaRKQuEsT was de- 
signed to predict the comoving number density of haloes at a given 
halo mass; i.e., the cumulative halo mass function, n(> Moo9,), for 
haloes with masses Mog9p > 10!2 h-!Mo. Also, Nishimichi et al. 


(6) 


7 We refer the reader to the original publications about the emulators for 
details about the adopted cosmological models and simulations. 
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Figure 10. Boost factor, B(k), comparison between the measurements from the mean of the MTNG740-DM-1-A/B pair of simulations (black lines) and the 
predictions from the CosmicEmu (orange lines), Bacco (blue lines) and EucLipDEMuLaTor2 (green lines) emulators at redshifts z = 0 (left panel) and z = 1 
(right panel). The lower subpanels show the relative difference between our simulations and the emulator predictions, while the dashed and dotted horizontal 


lines indicate a 1% and 5% difference, respectively. 
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Figure 11. Halo mass function comparison between the MTNG740-DM simulations and the CosmicEmtu (left panel) and the Dark Quest (right panel) emulators 
at z = 0 and z = 1, as labelled. The lower subpanels show the relative difference between our simulations and the emulator predictions, the dashed and dotted 


horizontal lines indicate a 5% and 10% difference, respectively. 


(2019) applied a correction to the halo mass given by, 
M200m = (1 — Npo??) M360 (7) 


where Np is the number of particles of a given halo, and a 
is the spherical overdensity mass definition obtained from the halo 
catalogue. We therefore apply the same correction to our M2o9m 
haloes to have a fair comparison. 

The measured HMFs from our halo catalogues at z = 0 and 
z = 1 are compared against the DarKQueEst predictions in the 
right panel of Fig. 11. In this case, we find a disagreement of 
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10% for haloes less massive than Mogop = 10!3-5 h-'Mo. How- 
ever, the differences are within the 5% level in the mass range 
10!3-5 < Mogop/[h7!Mo] < 10!4-5. The large deviations at the 
low-mass end are potentially due to the different halo finder used by 
Nishimichi et al. (2019) to construct their halo catalogues. DARK- 
Quest was built and validated with Rockstar catalogues (Behroozi 
et al. 2013), while our MTNG haloes were identified by the FoF algo- 
rithm combined with SuBFinp. Previously, Knebe et al. (2011) have 
shown that differences in the abundance of haloes of up to ~ 10% 
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Figure 12. The ratio between the measured non-linear power spectrum of the MTNG740-DM simulations and the smooth (no-wiggle) linear theory power 
spectrum from z = 30 (upper left panel) to z = 0 (lower right panel). Note that we are showing the mean over the A and B runs (red line), while the green line 
in each panel shows the ratio between the linear theory P(k) with respect to its smooth counterpart for each selected redshift. The vertical dashed line indicates 
the fundamental mode of the box, kpox, given by Eq. (2). In addition, we show the non-linear prediction from Havorir at redshifts z = 1 and z = 0. 


can occur between Rockstar and SUBFIND, which is consistent with 
our findings in the lower-right panel of Fig. 11. 


5 LARGE-SCALE MATTER AND HALO CLUSTERING 


One of the main goals of galaxy surveys is to measure the cluster- 
ing of matter and galaxies on the largest scales, r > 150 h-'Mpc 
(k<O.1h Mpc7!), to accurately determine the position of the bary- 
onic acoustic oscillations (BAO) peaks with an unprecedented level 
of precision. In this Section, we present an analysis of the large-scale 
clustering of matter and dark matter haloes. Note that we are us- 
ing the average of the MTNG740-DM-A and -B realisations, unless 
otherwise stated. 


5.1 Matter and halo clustering at the BAO scale 


The volume of the MTNG740 simulations allows us to compute 
the dark matter power spectrum with high accuracy at the BAO 
scales down to the present time. Fig. 12 displays the measured matter 
power spectrum from our MTNG740-DM simulations divided by the 
smoothed (no-wiggle) linear theory prediction from z = 30 (top left 
panel) to z = 0 (bottom right panel). At very early times (z = 30), the 
measured power spectrum matches perfectly with the linear theory, 
where we can identify up to six BAO peaks. We start to appreciate 
small deviations with respect to the linear predictions on small scales, 
k > 0.3hMpc™!, at z = 15 of order of 1 percent. At intermediate 
redshifts, z = 7 and z = 3 (middle row of Fig. 12), the non-linear 


effects are visible beyond the position of the third (k ~ 0.2 h~'!Mpc) 
and second (k ~ 0.13 h Mpc!) BAO peak, respectively, showing an 
enhancement of 5% and 15% in the clustering signal at the smallest 
scales. At low redshifts (z < 1; see the bottom row of Fig. 12), 
we detect negative and positive shifts of the BAO peaks of order 
2 — 3% due to non-linear effects, extending to the largest scales. 
For comparison, we show the predictions from the HaLorit model 
(Takahashi et al. 2012) at z = 1 and z = 0 (see the blue solid lines in 
the lower panels of Fig. 12). We find that HALorir predicts the same 
shifts of the BAO peaks as found in the simulation measurements. 

We now analyse the large-scale halo clustering in Fourier space. To 
do so, we select haloes according to their mass (M299,) and velocity 
(Vmax); note that we only use central subhaloes (selected as the largest 
bound structure identified by SuBFrinp in every FoF halo) in our Vinax 
samples. The maximum of the circular velocity, Vmax, can be used 
as a proxy to identify galaxies in dark matter subhaloes and make 
mock catalogues through the subhalo abundance matching technique 
(SHAM; see e.g., Conroy et al. 2006; Gerke et al. 2013; Klypin 
et al. 2013; Reddick et al. 2013). In each case, we select two halo 
samples with different space density, ny = 1 x 10-3 h>Mpe~? and 
ny = 6.86 x 10-4 h>Mpc™, at z = 1. We select samples with these 
specific number densities at a redshift of unity because this makes 
them close to the expected spatial density of DESI OII emission line 
galaxies (ELGs) and Euclid Ha emitters (DESI Collaboration 2016; 
Blanchard et al. 2020), respectively. 

We measure the monopole of the redshift-space halo power spec- 
trum because we want to draw a more direct analogy between the halo 
and galaxy clustering. Also, studying the clustering of Vinax-selected 
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Figure 13. The measured redshift-space dark matter halo power spectrum at z = | (blue dots with error bars) for two different M2o9-— (upper panels) and 
Vinax Selected (lower panels) halo samples with number densities nj, = 1 x 1073 hA>Mpc™ (left panels) and ny, = 6.86 x 107+ h3Mpce™? (right panels), as well 
as the best-fitting model (black solid line). The error bars correspond to the standard deviation over six Po(k) measurements. 


haloes is similar to studying the clustering of mock galaxies from a 
SHAM catalogue. We here use the distant-observer approximation to 
shift the positions of haloes from real- to redshift-space, we use the 
three coordinates, %, $ and Z, as the line-of-sight (LOS) to generate 
three redshift-space catalogues for one real-space catalogue where 
the new coordinates are, 


(1+z)vq 
H(z) e\|> 


where r is the comoving coordinate vector in real space, s is the 
equivalent of this in redshift-space, and z is the redshift. H(z) is 
the Hubble parameter, v and é) are the components of the peculiar 
velocity and the unit vector along the LOS direction. So, in total 
we have six redshift space catalogues. Then, the monopole of the 
redshift-space power spectrum can be obtained by 


s=rt+ 


(8) 


1 
Po(k) = I Pk, w) dye, (9) 
where P(k, j2) is the full two-dimensional power spectrum, and sz is 
the cosine of the angle between the separation vector, s or k, and the 
line-of-sight in configuration or Fourier space. 

We use NBopykirt to measure the halo power spectrum from the 
simulation for wavenumbers ky, < k/[hMpc!] < kny, where 
Apox is the fundamental mode of the box given by Eq. (2) and kny is 
the Nyquist frequency computed by Eq. (5) with Nerid = 256 (which 
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is here the mesh resolution of the power spectrum measurement), 
using linear k-bins with width Ak = 0.005 hMpc~! and 30 linearly 
spaced bins between 0 and 1 for w. 

The halo clustering measured from the MTNG740-DM simula- 
tions is shown in Fig. 13 (blue symbols with error bars), for each halo 
sample at z = 1. We focus only on large scales, k < 0.5 hMpc"!, 
where we can clearly see the BAO features in the power spectrum. 
Also, we can see that samples with a lower number of objects; i.e., 
np = 6.86 x 10-4 h>Mpc->, have a higher clustering amplitude than 
the larger samples (ny = 1 x 1073 h?Mpe~3). This is because the 
low-density sample is influenced more strongly by massive objects 
which are more highly clustered. 

In cosmological analyses of the BAOs, it is common to extract the 
peak position through the dilation parameter, a, which is related to 
the spherical average distance (Eisenstein et al. 2005). For a > 1 the 
peak is moved to smaller scales, while a < 1 shifts the peak to larger 
scales (Angulo et al. 2008; Anderson et al. 2014; Ross et al. 2015). 

As a proof-of-concept example, we estimate the dilation parameter 
from our halo samples by fitting the monopole of the power spectrum 
from our simulations to a simple template used by Springel et al. 


(2018), 
PUY (k) = BY, (A + Agk + A3k7)Phin(k/a) . (10) 


where Pji, is the linear theory power spectrum, Bp is a bias pa- 
rameter, and A,, Az and A3 are further free parameters. We fit our 
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Figure 14. Scale dependent halo bias, b(k) = V Py(k)/Pm(k), for Moo9-—selected and Vinax—selected haloes at z = 0 (left panel) and z = 1 (left panel), for 
the MTNG740-DM-| simulations. Different line-styles represent different types of simulations, as labelled. 


measurements on scales k < 0.5hMpc™!. The best fit model for 
each sample is shown together with the measurements in Fig. 13. In 
all cases, we can constrain the dilation parameter with a precision of 
Ca/a ~ 0.08 — 1.4 percent. 


5.2 Linear halo bias 


Haloes are biased tracers of the dark matter density field, hence 
the relation between the distribution of haloes and matter can be 
described by the linear halo bias b, defined as 


b = oy/6, (11) 


where 6}, is the halo density contrast and 6 is the density contrast 
of matter. In terms of a Fourier-space analysis of the clustering, the 


halo bias can be estimated by, 


Py(k) 
Pm(k)’ 


b(k) = (12) 
where P,(k) and Pm(k) are, respectively, the halo and matter power 
spectrum in real space. At sufficiently large (linear) scales, the halo 
bias is expected to be a scale independent constant, and its amplitude 
should depend mostly on halo properties, such as the halo mass. 
This linear halo bias is a basic ingredient of the halo model (see 
e.g., Cooray & Sheth 2002). Since haloes host different galaxy types, 
it is important to accurately describe the halo bias in order to use 
it to model galaxy clustering and the galaxy-halo connection (for a 
review, see Wechsler & Tinker 2018). 

We measure the halo bias for fixed My99¢ and Vmax halo samples 
at z = O and z = 1, as shown in Fig. 14. We work with haloes with 
masses above Myg9. = 10!! A-'Mo (M1), 10!! h-'Mo (M2), 
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Figure 15. Matter power spectra for different mass components of the MTNG740 hydrodynamic simulation (dashed lines) at redshifts z = 0 (right panel) and 
z = 1 (left panel). We include the measurement of the TNG300 (dotted lines) and TNG100 (solid lines) simulations, for comparison. The lower subpanels 
display the relative difference between the TNG300/TNG100 and MTNG740 measurements. 
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Figure 16. Impact of baryonic physics on the total matter spectrum at z = 0 
measured from the MTNG740 runs (red solid line) in comparison with the 
TNG100 (cyan dashed line) and TNG300 (brown dashed line) simulations. 
The dashed horizontal lines and the grey shaded region indicate a relative 
difference of 2 and 10 per cent, respectively. 


10!2 h-1Mo (M3), 10!25 A-!Mo (M4) and 10!3 A-!Mo (M5), as 
well as with samples with velocities higher than Vinax = 10!-8kms7! 
(V1), 10?kms7! (V2), 107-2 kms! (V3), 102-+kms7! (V4) and 
102-6 km s~! (V5). Similarly to the redshift-space power spectrum 
measurements discussed in Sec. 5.1, we measure the real-space clus- 
tering in Fourier space for all halo samples in the range kpox < 
k/[hMpc7!] < 1 with k-bins with width Ak = 0.005 h Mpc7!. We 
additionally show in Fig. 14 the individual measurements from the 
A (dashed lines) and B (dot-dashed lines) realisations, besides the 
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1.05 
— MTNG740 


Figure 17. Suppression of the matter power spectrum due to baryonic physics 
in the MTNG740 hydrodynamic simulation relative to MTNG740-DM-1 
(solid lines) and the predictions from the Bacco emulator (symbols) at z = 0 
and z = 1. The lower subpanel shows the relative difference between the 
simulation and the emulator. 


mean of the paired runs (solid lines). In general, more massive haloes 
and haloes with higher Vinax velocities display a larger halo bias. 
From the upper panels of Fig. 14 we can see that the bias of 
Mp 0¢-selected haloes shows an exclusion effect on small scales for 
the most massive samples (M4 and M5), where the bias drops at 
k~O0.5h Mpc™! at both output times, z = 0 and z = 1. Also, the bias 
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Figure 18. The differential halo mass functions as a function of halo mass (M299-; left panel), and the Vinax subhalo functions (right panel) measured from the 
MTNG740 hydro and the MTNG740-DM-1 simulations at z = 0, as labelled. The lower subpanels show the ratio between the hydro halo mass and velocity 
functions with respect to their dark matter-only versions. We also show the mean of the A and B realisations as a red line. 


at early times, e.g., z = 1 (upper right panel), has a higher amplitude 
for fixed halo mass samples than their present time counterparts 
(upper left panel). In all cases, we find that the halo bias is constant 
on very large scales, k < 0.1 hMpc"!. However, on mildly non- 
linear scales, scale-dependent effects set in and become appreciable 
on scales k > 0.3 hMpc™!. 

On the other hand, velocity selected halo samples (bottom panels 
of Fig. 14) show qualitatively similar results on large-scales than 
Mo 9¢-selected haloes. However, we do not observe exclusion ef- 
fects for samples with higher Vmax, i.e., V4 and V5. But we observe a 
strong scale-dependent bias on scales k > 0.3h Mpc™!. Our simula- 
tions make very accurate quantitative predictions for these non-linear 
effects in dark matter-only models, but we caution that they can be 
affected by baryonic physics impacts, which constitute the primary 
systematic uncertainty. 

Note that the A- and B- runs produce similar results, but the 
average of both realisations tends to cancel out fluctuations on the 
largest scales, producing a much smoother measurement of the linear 
halo bias. 


6 BARYONIC EFFECTS ON BASIC MATTER AND HALO 
STATISTICS 


So far we have focused our attention on the results from the MTNG 
DM-only runs, ignoring baryons, as well as neutrinos. In this section, 
we briefly discuss the impact of baryons on basic matter and halo 
statistics. Previously, Springel et al. (2018) presented a high-level 
study of the matter and halo clustering results of the full-physics 
IllustrisTNG simulations, giving some hints of the effects of baryonic 
feedback of the TNG model on basic matter and halo statistics. In a 
companion paper, Pakmor et al. (2022) explores the consistency of 
the MTNG740 full-hydro run with the smaller box TNG simulations. 
Here we extend this analysis and consider the baryonic impact on 
basic halo and matter clustering statistics, while a more general study 


of the baryonic effects also on other observables, such as weak lensing 
and general galaxy statistics, will be presented in forthcoming work. 


6.1 Clustering of matter components 


In Figure 15 we show the clustering power spectra of different matter 
components in our large MTNG740 hydro run at redshifts z = 0 and 
z = 1; i.e. separately for dark matter, gas, stars, and black holes. We 
also compare to equivalent measurements for TNG100 and TNG300 
from the IllustrisTNG project (Springel et al. 2018). We find in gen- 
eral good consistency between the MTNG740 and the TNG100/300 
simulations, with some small differences that can be understood 
as arising from differences in mass resolution, and, in particular, 
in simulated volume. Reassuringly, there is excellent agreement on 
large scales, with MTNG740 extending the baryonic results to much 
larger spatial scales, throughout the BAO region, than possible for 
IlustrisTNG. 

As shown in Fig. 16 (and in our companion paper by Pakmor et al. 
(2022, their Fig. 4)), the suppression of the total matter clustering 
predicted by MTNG740 is also in nearly perfect agreement with that 
found in TNG100, confirming that there is a reduction of power of 
around 20% at k ~ 10h Mpc™! due to AGN feedback, and a strong 
enhancement of the clustering on small scales (k > 50hMpc~!) 
compared to dark matter-only models as result of the dissipative 
formation of galaxies. We note that the MTNG740 suppression is 
almost identical to the TNG300 prediction at k < 5h Mpc™!. 

This suppression needs to be taken into account, in particular, 
in weak gravitational lensing analysis, unless one restricts them to 
comparatively large scales and thus gives up on considerable cosmo- 
logical information. It is thus important to develop tools to accurately 
model the suppression. In Figure 17, we consider how well this sup- 
pression can be described by the Bacco emulator (Arico et al. 2021) 
at z = O and z = 1. The emulator matches our numerical results 
with impressive accuracy, albeit limited to the large-scale regime, 
k<5 hMpc"!, for which the emulator has been calibrated. Note 
that we use the Bacco parameters that match with the TNG300 sup- 


MNRAS 000, 1—23 (2022) 


20 ~~ C. Herndndez-Aguayo et al. 


z=0 


Moo0c — selected haloes 


MTNG740 
MTNG740-DM 

np, =3x 10 3n3Mpc 3 
Np =1x 10-3h?Mpe~* 
Np =3 x 10-*h?Mpe~ 


Aé(r)/&(r) pmo 


10° 10! 


£o(s) 


MTNG740 
MTNG740-DM 

np =3 x 10°A3 Mpc? 
nyp=1x 10-°h3Mpe? 
ny =3 x 10 *h?Mpe? 


A&o(s)/€0(s)pmo 


s [h~*Mpc| 


MTNG740 
MTNG740-DM 
10? np =3 x 10-3? Mpe3 
np =1x 10-A?Mpe? 
np, =3 x 10-*h3Mpe3 


z=0 


Vinax — Selected haloes 


10? 
z=0 


Vinax — Selected haloes 


10? 


MTNG740 
MTNG740-DM 

np =3 x 10h? Mpc? 
nyp=1x 10-°A® Mpc? 
ny, =3 x10 *h?Mpe ® 


107 


eoofor 
CNBRDAWO 


A&o(s)/€0(s)pmMo 


s [h~ Mpc] 


Figure 19. The halo two-point correlation function in real-space (upper panels) and redshift-space (lower panels) for three different halo number density cuts, 
npn = 3x 10-3 h3>Mpc™ (blue lines), ny = 1 x 10-3 h3Mpc~3 (green lines) and np, = 3 x 10+ h>Mpc-3 (red lines), at redshift z = 0 for Mzg9-— and 
Vinax —Selected haloes, as labelled. The lower subpanels show the relative difference between the MTNG740 hydro and the MTNG740-DM measurements. 


pression obtained by Arico et al. (2021), which seems adequate as 
MTNG740 and TNG300 give almost identical results (Pakmor et al. 
2022) in the k-range of the Bacco emulator. The level of agree- 
ment is nevertheless encouraging and suggests that parametrisations 
like Bacco can be of great help to marginalise over complicated 
astrophysical effects such as AGN feedback. 


6.2 Halo abundance and clustering 


Previously, it has already been shown that the presence and evolution 
of gas and stars in hydrodynamical simulations affects the mass 
of haloes (Vogelsberger et al. 2014; Springel et al. 2018). This is 
because gas can be expelled from the inner regions of haloes thanks 
to the baryonic feedback processes, directly modifying the mass of a 
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halo and affecting its ability to grow through the attraction of further 
matter. This mass change is expected to not only affect the abundance 
of dark matter haloes as a function of mass, but also their clustering. 


In the left panel of Fig. 18, we show the differential halo mass func- 
tion (HMF) measured from the full-hydro (solid lines) and DM-only 
(dashed lines) runs of the MTNG740 (red lines), TNG100 (green 
lines) and TNG300 (orange lines) simulations at z = 0, for com- 
parison. In the bottom panels we show the residuals in both cases. 
Compared to dark matter only, we can see that the MTNG740 hydro 
run predicts 20% fewer low-mass haloes (M290¢ < 10!! 7-!Mo), an 
additional 10% deficit of haloes with masses Mzgo¢ ~ 10!3 h7~!Mo, 
and also an enhancement of the most massive haloes (Moo9- > 
10!4 2-!Mo) of ~ 5%. These results reaffirm the findings of Springel 
et al. (2018) for the TNG physics models. 
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Full hydrodynamical simulations also offer the possibility to ex- 
plore predictions for the number of dark matter substructures. In the 
right panel of Fig. 18, we show the differential subhalo Vinax func- 
tion at z = 0. We use only the most massive (central) subhaloes to 
measure such velocity functions from the full-physics and DM-only 
MTNG740 runs, as well as for TNG100 and TNG300. 

We find that there is a lack of low-velocity subhaloes, Vmax ~ 
125km sl, of around 20% in the hydro simulation; this is consistent 
with what we showed in the right panel of Fig. 18 for low-mass 
haloes. Also, we can see that there is a peak in the MTNG740- 
hydro velocity functions at Vmax ~ 310km s! producing 30% more 
subhaloes with respect to DM-only simulations. This peak divides 
the subhalo population into two samples. On one hand, we have low- 
velocity subhaloes that are affected strongly by supernovae feedback 
and which are related to low-mass haloes, and on the other hand, 
there are high-velocity subhaloes that occupy cluster-like haloes, and 
which are more strongly influenced by AGN feedback. 

Of particular relevance observationally is how the structural 
changes in the properties of haloes are reflected in their clustering 
properties. Note that changes of halo properties invariably intro- 
duce differences in sample selection, which can in turn influence 
the clustering signal. To obtain a first quantitative assessment of this 
we consider in Figure 19 halo samples of different space densities 
(np = 3x10~3 h3>Mpe~3, 1 x1073 A2Mpe~? and 3 x10~4 A3Mpe~), 
either selected as a function of halo mass, or as a function of the max- 
imum circular velocity. Note that the latter can be strongly affected 
by effects such as adiabatic compression due to baryonic physics, 
even if the halo virial mass itself does not change (as we showed in 
the left panel of Fig. 18). We give measurements at z = 0 both for 
real-space and redshift-space. 

Interestingly, we find substantial differences in the clustering am- 
plitudes even on large scales. While they are comparatively small at 
the level of 2 percent for mass selected samples in real-space, they 
become larger and sample-density dependent in redshift space. For 
circular velocity selected samples, the differences are substantial, 
amounting to 50% at large scales, approximately independent of the 
sample density, and whether real-space or redshift-space is consid- 
ered. For distances of order 1 h-'!Mpc and smaller, halo exclusion 
effects that depend on the sample density become clearly notice- 
able. Overall, these results stress once more the strong impact of the 
sample selection on any clustering analysis (see also Angulo et al. 
2012). 


7 SUMMARY AND CONCLUSIONS 


In this paper, we have introduced the MILLENNIUMTNG simula- 
tion project. Our flagship runs consist of two fixed-and-paired DM- 
only simulations with 43203 dark matter particles in a volume of 
500° n-3 Mpc? (which is the same volume as in the iconic MILLEN- 
NIUM run; Springel et al. 2005), a single full-physics hydrodynamical 
simulation which evolves 2 x 4320% dark matter and gaseous reso- 
lution elements in the same cosmological volume from z = 63 to 
the present time, and an extremely large simulation with more than a 
trillion dark matter particles which also includes massive neutrinos 
(Xm, = 0.1 eV) in a volume of (3 Gpc)?. The full-hydro simulation 
employs the ILLustrisTNG model of galaxy formation (Weinberger 
etal. 2017; Pillepich et al. 201 8a). In addition, we have complemented 
the simulation set with a number of smaller runs to investigate res- 
olution dependence, box-size dependence, as well as the impact of 
different neutrino masses. 

In this paper, we have introduced the primary technical aspects of 


the runs (see Table 1) and describe their data products. However, our 
analysis only focuses on the matter clustering and halo statistics of the 
DM-only and full-physics MTNG740 runs. An in-depth analysis of 
the MTNG-neutrino simulations will be presented in a forthcoming 
paper. 

We have performed a series of convergence tests by comparing 
the measured non-linear matter power spectrum (see Fig. 2), the 
halo mass functions for two halo mass definitions (M99. and Mfor; 
shown in Figs. 4 and 5, respectively), and the real-space halo clus- 
tering for three halo mass selected samples (Fig. 6) at z = O and 
z = 1, of our level-1 run with the lower resolution simulations. We 
find sub percent agreement between the 4320? and 2160? resolutions 
for all measurements over a large range of scales and halo masses. 
Additionally, we have assessed the convergence of our numerical 
calculations by comparing the non-linear matter clustering and halo 
abundances of the MTNG740-DM and MTNG185-DM boxes (see 
Fig. 7). The two simulated volumes agree very well (< 2%) and 
the largest differences are directly associated with effects due to the 
limited volume of the MTNG185-DM run. 

We have shown that computing two realizations of a given model 
using fixed-and-paired N-body simulations allows us to make accu- 
rate cosmological predictions at large scales despite the still quite 
limited volume (see Fig. 8), and that this technique yields results 
comparable to those predicted by the mean of one hundred inde- 
pendent realisations (Fig. 9). This confirms that the fixed-and-paired 
technique effectively boosts the statistical power that is achievable 
with a given simulation volume. 

We have compared the non-linear boost factor, Eq. (6), and the 
halo mass functions from a set of recent cosmological emulators 
(CosmicEmu, Bacco, EUCLIDEMULATOR?, and DAaRKQUEST) to mea- 
surements from our simulations, finding reasonably good agreement 
in all the cases over the range where the emulators were calibrated 
(Figs. 10 and 11). This can be viewed both as a reassuring confir- 
mation of the high-quality calibration reached by these emulators, as 
well as an independent test of our simulation results if one already 
trusts these emulators. 

We have also given a first impression of the large-scale clustering 
of matter and haloes around the BAO scale in Fourier space (Figs. 12 
and 13). We have shown that the MTNG740 simulations have enough 
statistical power and volume to extract the BAO features through the 
dilation parameter, a. We also show the linear halo bias for Maoo¢ 
and Vmax-selected haloes in Fig. 14. 

We have validated the clustering of different matter components 
(dark matter, gas, stars and black holes) of the full-physics MTNG740 
simulation by comparing with results from the TNG100 and TNG300 
simulations (see Figs. 15 and 16), finding good agreement between 
MTNG and its TNG predecessors. We also compared the suppression 
of the total matter power spectrum due to baryonic physics between 
MTNG and the Bacco emulator (see Fig. 17), finding sub percent 
agreement on scales k < 5hMpc7!. The halo and substructure 
abundances (Fig. 18) are also in good agreement between MTNG740 
and TNG100/300. 

Finally, we have explored baryonic effects on halo clustering in 
real- and redshift-space at z = 0 in Fig. 19. We find that baryons 
affect the clustering of Mz99--selected haloes between ~ 2% in real- 
space up to > 20% in redshift-space. The clustering of Vmax-selected 
haloes is modified even more, here the inclusion of baryons causes 
differences of ~ 50% in both real- and redshift-space. 

Numerical calculations of the size of the MTNG simulations are 
needed to assist modern galaxy surveys in achieving their goals to 
constrain cosmological parameters at the sub percent level. With the 
MING project, we aim to produce realistic mock galaxy catalogues 
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using the physically motivated TNG model in combination with semi- 
analytical models of galaxy formation. However, before presenting 
galaxy mock catalogues constructed in this manner, we need to care- 
fully assess the reliability of the model foundations; i.e. whether the 
MTNG simulations pass very high standards of numerical accuracy. 
In this paper, we have examined this question and arrived at an affir- 
mative assessment, making us confident that the MTNG simulations 
are a powerful tool to study a number of cosmological questions, 
some of them we address both in our companion papers and in our 
forthcoming work. 
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